NIST-GCR-93-637 


A COMPUTATIONAL MODEL FOR THE 
RISE AND DISPERSION OF WIND- 
BLOWN, BUOYANCY-DRIVEN PLUMES. 
PART Il. LINEARLY STRATIFIED 
ATMOSPHERE 


Xiaoming Zhang and Ahmed E Ghoniem 


Massachusetts Institute of Technology 
Cambridge, MA 02139 


NUT 


United States Department of Commerce 
Technology Administration 
National Institute of Standards and Technology 


NIST-GCR-93-637 


A COMPUTATIONAL MODEL FOR THE 
RISE AND DISPERSION OF WIND- 
BLOWN, BUOYANCY-DRIVEN PLUMES. 
PART II. LINEARLY STRATIFIED 
ATMOSPHERE 


Xiaoming Zhang and Ahmed E Ghoniem 
Massachusetts Institute of Technology 
Cambridge, MA 02139 


December 1993 


Sponsored by: 

U.S. Department of Commerce 

Ronald H. Brown, Secretary 

Technology Administration 

Mary L. Good, Under Secretary for Technology 
National Institute of Standards and Technology 
Arati Prabhakar, Director 


Notice 


This report was prepared for the Building and Fire 
Research Laboratory of the National Institute of 
Standards and Technology under grant number 60NANBOD1036. 
The statements and conclusions contained in this report 
are those of the authors and do not necessarily reflect 
the views of the National Institute of Standards and 
Technology or the Building and Fire Research Laboratory. 


til 


Prepared for publication at Atmospheric Environment 
July, 1993 


A COMPUTATIONAL MODEL FOR THE RISE AND DISPERSION OF WIND-BLOWN, 
BUOYANCY-DRIVEN PLUMES 
PART I. LINEARLY STRATIFIED ATMOSPHERE 


Xiaoming Zhang and Ahmed F. Ghoniem 
Department of Mechanical Engineering 
Massachusetts Institute of Technology 

Cambridge, MA 02139 


Corresponding Author: 


Ahmed F. Ghoniem, Professor 
Massachusetts Institute of Technology 
77 Massachusetts Avenue, Room 3-342 
Cambridge, MA 02139 

Phone: (617) 253-2295 

FAX: (617) 253-5981 

email address:ghoniem@mit.edu 


Key words; buoyant plumes, vortex simulation, stratified atmosphere. 


A COMPUTATIONAL MODEL FOR THE RISE AND DISPERSION OF WIND-BLOWN, 
BUOYANCY-DRIVEN PLUMES 
PART I. LINEARLY STRATIFIED ATMOSPHERE 


Xiaoming Zhang and Ahmed F. Ghoniem 
Department of Mechanical Engineering 
Massachusetts Institute of Technology 

Cambridge, MA 02139 


Abstract 


A multi-dimensional computational model of wind-blown, buoyancy-driven flows is 
applied to study the effect of atmospheric stratification on the rise and dispersion of plumes. The 
model utilizes Lagrangian transport elements, distributed in the plane of the plume cross section 
normal to the wind direction, to capture the evolution of the vorticity and density field, and 
another set of elements to model the dynamics in the atmosphere surrounding the plume. 
Solutions are obtained for a case in which atmospheric density changes linearly with height. 
Computational results show that, similar to the case of a neutrally stratified atmosphere, the 
plume acquires a kidney-shaped cross section which persists for a long distance downstream the 
source and may bifurcate into separate and distinct lumps. Baroclinic vorticity generated both 
along the plume boundary and in the surroundings are used to explain the origin of the distortion 
experienced by the plume and inhibiting effect of a stratified atmosphere, respectively. The 
vorticity within the plume cross section forms two large-scale coherent eddies which are 
responsible for the plume motion and the entrainment. Prior to reaching the equilibrium height, 
the computed plume trajectory is found to follow the two-thirds law, when extended to include 
the initial plume size, reasonably well. The entrainment and the added mass coefficients, 0.49 
and 0.7, respectively, are obtained from the numerical results over a wide range of the buoyancy 
ratio, defined as the ration between the plume buoyancy and the degree of background 
stratification. In the case of strong stratification, the plume trajectory shows weak, fast decaying 
oscillations around the equilibrium height. The origin and decay of these oscillations are 
explained using a simple analytical model. 


Nomenclature 


B Buoyancy Ratio; 
4 Gravitational acceleration; 
h Spatial discretization length in the plume-air interface; 
H Vertical dimension of the computational domain; 
hs Spatial discretization length in the stratified background; 
HT Initial plume height above the ground; 
ky Added mass coefficient; 
N Brunt-Vaisaila buoyant frequency; 
N’ Modified buoyant frequency; 
R Square root of the plume cross sectional area; 
Ro Initial Radius of the circular plume cross section; 
Req Equivalent radius of plume cross section; 
Ry Major (horizontal) axis of elliptical plume cross section; 
Rz Minor (vertical) axis of elliptical plume cross section; 
Rep Buoyancy Reynolds number; 
t Same as the normalized x coordinate; 
U Homogeneous wind speed; 
V = R g , plume buoyancy velocity; 
w Plume vertical velocity; 
x Horizontal wind direction; 
y Horizontal direction normal to the wind direction; 
z Vertical direction, plume rise; 
B Entrainment coefficient; 
Pa Arr density; 
Po Reference density; 
Pp Plume density apart from reference density; 
Ps Background air density apart from reference density; 
E Vertical displacement of plume from its equilibrium height; 
At Time step; 
E = eet plume mass flux ratio; 
(2) 
Superscript 


* 


Dimensional quantities. 


1. Introduction 

The dispersion of strongly buoyant plumes, similar to those generated from massive fires, 
cooling towers and tall industrial and power plant stacks, in a stable atmosphere is important in 
many practical applications. Stratification can limit the plume rise and distort its dispersion 
patterns in ways that worsen its impact on local air quality. Current models (Briggs, 1975) used 
to describe plume rise and dispersion in a stratified atmosphere, are based mostly on assumptions 
regarding entrainment and self-similarity and involve several adjustable constants (see Section 3 
for a review). Application of computational methods is impeded by the difficulties encountered 
in the modeling of the buoyancy-generated turbulence, the large numerical diffusion associated 
with Eulerian discretization schemes, and the computational cost required when dealing with 
large domains. 

A recently developed computational model, based on the Lagrangian interpretation of the 
dynamics of buoyancy-driven flows, attempts to overcome some of these difficulties (Zhang and 
Ghoniem, 1993, hereafter referred as Part I). The numerical model utilizes the vortex element 
and transport element methods to cove the equations governing a wind-blown, buoyancy-driven 
plume derived directly from Navier-Stokes equations. This paper describes an extension of this 
model, originally intended for flows in a neutrally stratified atmosphere, to incorporate the effect 
of atmospheric density stratification on buoyant plume motion. We focus on the case in which 
the potential density! (referred to simply as density) of the air decreases linearly upwards. This 
form of stratification occurs frequently in elevated layers during the day and in the lowest 100 m 
or so at night (Briggs, 1975). In a following publication, the effect of an inversion iayer will be 
described in detail. 

Atmospheric turbulence may be neglected if one is interested in buoyant plumes not too 
far away from the emission source, where the plume buoyancy dominates the motion. In general, 


and depending on how it is generated, there are two types of atmospheric turbulence. The first, 


1 Potential density is defined as the density an air parcel would reach if it was brought adiabatically to a 
standard pressure, which for convenience is usually taken as the pressure at the ground. A neutrally 
stratified atmosphere is one in which the potential density is constant. 
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characterized by Vu”; where u’ = u - U is the fluctuating component of the instantaneous wind 


velocity, u, around its mean value, U, is due to wind shear instability. On the other hand, the 


plume buoyancy can be characterized by a buoyant velocity, V = , where Po, Pp, RF, 
oO 


and g are the reference density, the plume deficient density, the characteristic length of the plume 
cross section, and the gravitational acceleration, respectively (see Part I). A criterion for 


neglecting this type of atmospheric turbulence may be given by 


(1) 


The second type of atmospheric turbulence is associated with the convective instability 
due to the uneven heating of air parcels in the atmosphere. A velocity scale representing the 


intensity of this turbulence over a horizontal distance equal to the plume size is given by 


, where Af, is the scale of air density variation over a horizontal distance R. A 
Oo 


corresponding criterion for neglecting this convective atmospheric turbulence is thus 


in >>ak. (2) 
Apa 


In this paper, we neglect atmospheric turbulence and focus on plume rise and dispersion 
due to buoyancy and buoyancy-generated turbulence in a stratified atmosphere (the effect of the 
former will be incorporated in future efforts). In Section 2, we describe the physical mechanism 
of baroclinic-vorticity generation in buoyant flows, and summarize how atmospheric 
stratification is implemented in our model. In Section 3, we review the integral models which 
have been developed to characterize the plume trajectory, and refer to field observations on 
plume dispersion in a linearly stratified atmosphere. In Section 4, we present the computational 
results in the form of the plume trajectory and dispersion pattern, and compare our results with 
experimental data. An algebraic formula for the plume rise is obtained by modifying the 
classical two-thirds power law. The adjustable empirical constants in the formula: the 


entrainment coefficient, and the added mass coefficient are calibrated using the numerical results. 


y 


The paper ends with a brief conclusion in Section 5, and a note on plume motion near its 


equilibrium height in the Appendix. 


2. Background 
2.1. Vorticity Generation in a Disturbed Stratified Atmosphere 


Dynamically, plume rise and dispersion in a weakly turbulent atmosphere is driven by the 
potential energy of the plume material. Kinematically, the plume motion is driven by the 
vorticity generated due to the density difference between the plume material and the surrounding 
air (see Part I for detail). The latter view has the advantage of allowing one to describe the 
turbulence generated within the plume cross section as well. To simplify the following 
discussion, we use the Boussinesq approximation which, in our previous work, has been shown 
to be accurate when the relative density perturbation is less than 0.1. With this approximation, 
baroclinic vorticity is generated only where a finite horizontal density gradient exists, i.e. along 
the distorted plume-air interface which at high Reynolds number, remains confined to a small 
subset of the entire domain. We also assume, as in Part I, that the density distribution inside the 
initial plume cross section is uniform. The effects of density non-uniformity within the plume 
cross section, which is expected to be of secondary importance, will be investigated in the future. 

In the presence of atmospheric density stratification, as the plume rises, the initially 
horizontal density isolines in the background are distorted. A density field with non-zero 
horizontal gradient component is thus established, and baroclinic vorticity is generated in the 
atmosphere. This is illustrated schematically in figure 1 for a crosswind section of a rising plume 
in a continuously stratified, stable atmosphere. As shown in the figure, vorticity generated in the 
background induces a net downward motion on the plume in the direction opposite to the flow 
induced by the vorticity generated due to the plume's own buoyancy. Eventually, the plume 
reaches a terminal, equilibrium, or level-off height where the effect of the background vorticity 
balances that of the vorticity around its surface. Occasionally, the plume may experience weak, 
fast decaying oscillations around this height (Briggs, 1975) due to a mechanism that will be 


discussed later. It should be noted, as will be shown later, that even after the plume reaches its 


equilibrium height, vorticity generated both at the plume surface and in the surroundings still 


plays an important role in dispersing its material. 
2.2. Model Formulation and Numerical Implementation 


The partial differential equations governing the flow of a steady, three-dimensional, 
wind-blown, buoyancy-driven plume were presented in Part I. Incorporating the atmospheric 
density stratification in the formulation does not change the governing equations, except for the 
fact that vorticity can now be generated in the background and its induced field must be 
incorporated in the solution. 

In principle, in a continuously stratified atmosphere, we must discretize the entire space 
since non-zero horizontal density gradients could be established everywhere. However, the 
plume induced motion, and thus the distortion of the background, decays rapidly away from the 
plume center?. Moreover, the far field effect on the plume motion is generally negligible, see 
Section 4.1. Thus, to account for vorticity generated in the atmosphere noting that the air density 
isolines remain almost horizontal in the far field, we need only consider an area several plume 
characteristic length, R, around the plume center . Furthermore, in practice, the density gradient 
in the background is much smaller than that across the plume-air interface. Thus, we discretize 
the background density field into transport elements whose size is larger than the size of the 
elements used to discretize the plume-air interface. The discretization is schematically illustrated 
in figure 2. The choice of two different sizes of elements reduces the computational load while 


maintaining the numerical accuracy. 


2 Note that the distortion of the atmospheric density isolines, induced by the plume vorticity, is 
proportional to the velocity gradient, Au = - Ar, where Tis the circulation and r is the distance away from 


the vorticity center, i.e., the effect of vorticity induced distortion decays as the square of the distance away 
from the center of the vorticity. 
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3. Review 


Plume rise in a linearly stratified atmosphere has been analyzed in the literature. 
Denoting the reference density, taken here as the air density at the height of the plume source, by 


Po, We write the density distribution as 


Po + Py(ty,2) inside the plume 
p= { (3) 


Po + P,(ty.2) outside the plume 
where Pp and ps are the density perturbation of the plume and the air, respectively. The 


background stratification is characterized by the Brunt-Vaisaila frequency, N2 = - dbs ida As it 


rises, the plume entrains denser air until it reaches an equilibrium height where its average 
density is close to that of its surroundings. Meanwhile, the plume spreads out horizontally and 
vertically. The trajectory of a point source plume between the source and a point close to the 
level-off height has been described by the two-thirds power law (Briggs, 1969, 1975, Fay, 1973, 
Weil 1988). In a uniform wind with speed U, the trajectory of a buoyant plume generated from a 
finite size source can be obtained using the same procedure as in the derivation of conventional 
two-thirds power law after modifying the initial conditions to account for the finite initial plume 
size. The derivation follows Weil’s (1988) and is similar to that described in Part I of this study. 
The resulting expression of the plume rise z* is: 
was 3BFp N' x* oxi : N x* 
ze al tt cos NX*_) +R; | -Ri\ for Net <q (4) 
B \LQ+k)N2U U U 

where R5 is the radius of the initial plume cross section, ky is the “added” or “virtual” mass 
coefficient which accounts for the momentum of the ambient fluid displaced by the plume as the 
latter rises and of the small scale motion inside the plume; N' = N (1 + k,)!/2 is a modified 


buoyancy frequency which takes k, into account; B is an entrainment coefficient used to describe 


the growth of the plume radius due to entrainment3; and Fy = Ro " Ug pe is the buoyancy flux. 
Oo 


3 Weil (1988) used a modified entrainment rate related to the added mass coefficient in his derivation. 
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The constants in equation (4) have been obtained using laboratory experiments and field 
observations of point-source plumes. However, as pointed out by Turner (1986), this two-thirds 
power "law" fails near the level-off region because the entrainment and similarity assumptions 
used in its development become invalid. Similar reasons, coupled with the difficulty in defining 
the plume center and the size of its cross section in different experiments (List 1982, Gebhart et 
al, 1984), have been used to describe the origin of the following unresolved problems: 

(1) Good correlation between the two-thirds power law and different experiments are achieved 
by adjusting the entrainment coefficient for each case, i.e., the coefficients are not truly 
universal. From laboratory experiments, Slawson and Csanady (1967) found B = 0.43, while 
Hewett et al. (1971) got the best fit with B = 0.71. From field experiments, Hoult et al (1969) 
found that 0.3 < B<1.0, while Fay et al. (1970) obtained 0.5 < B < 2.0. 
(2) In models which assume that the plume cross section remains circular, the radius increases 
with the plume rise as 

R*=Ro + Bre. (5) 
Using data from the same experiment to fit the trajectory, equation (4), or the radius, equation 
(5), Briggs (1975) found B = 0.6 or 0.4, respectively. This inconsistency was also observed by 
other authors; Weil (1982) obtained B = 0.6 from the trajectory data while B =0.5 from the 
radius data. 

In the early development of the two-thirds power law, the added mass coefficient, ky , was 
not included (Briggs, 1969, Fay et al, 1970). Later k, was introduced (Briggs, 1975, Weil, 1988) 
to improve the accuracy of the trajectory law derived from the integral formulation. The added 
mass concept was originally developed for a solid body accelerating in a fluid (Yih, 1977) to 
account for the work done to increase the momentum of the ambient fluid. In the case of a 
circular cylinder moving in an inviscid flow, ky = J. The same value is usually suggested for 
plumes where the driving force is buoyancy. However, the actual value may deviate from unity 
since the plume is not a solid.body and its cross section is not circular. If fact, the plume cross 


section, as shown in Part I, is shaped in the form of an inverted kidney which evolves as the 
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plume rises. Furthermore, not only a significant portion of surrounding air is induced into 
motion due to the presence of the plume, but also motion at scales smaller than the plume radius 
is developed within its cross section. In practice, k, can be regarded as merely another empirical 
constant whose value should be adjusted to make the two-thirds power law fit the experimental 
trajectory better. 

The shape of the plume cross section is an important feature of plume dispersion, 
especially for large-scale, strongly buoyant plumes. As the plume rises, two counter-rotating 
vortices, which divide the plume material into two lumps, are generated. Plume bifurcation is a 
frequently observed yet not well understood phenomenon. Some of the laboratory and field 
experiments describing this bifurcation have been reviewed in Part I. Here we cite three more 
experiments, which have been documented carefully, in which bifurcation have been cited. The 
first is the fire plume experiment of Church et al (1980), conducted for the purpose of 
investigating the mechanism leading to the formation of large-scale vortices in fire plumes. 
These authors hypothesized that the generation of a double-vortex structure is due to the tilting 
and stretching of the vorticity originally associated with background wind shear, and the vorticity 
generated within the plume by the action of buoyancy and drag forces. The second is the Lidar 
measurements of cooling tower plumes (Hawley, 1985, Bennett et al, 1992). Sykes et al (1986) 
suggested that the interaction between the crosswind and the initial vertical momentum of the 
plume may be the cause for the bifurcation. The third experiment is the laboratory water plume 
model of Alton et al (1993). Clearly, the mechanism leading to plume bifurcation requires 


further investigation (Zhang and Ghoniem, 1993). 
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4. Numerical results 
4.1. Entrainment, added mass and plume trajectory 


The notations and normalization procedure used in this paper are the same as those 
employed in Part I. To characterize the background in which the plume disperses, we define the 


buoyancy ratio as the ratio between the plume buoyancy to the degree of stratification, 
fe A Se EE (6) 
N?R2 N?R 
where € = be is the plume deficient mass flux ratio and R is the plume length scale, taken as the 
oO 


square root of the initial plume cross sectional area. In terms of the buoyancy ratio, equation (4) 


becomes 
z a [sim Ate cos Tasers Lia i), (7) 


where the effects of plume buoyancy and background stratification have been grouped into a 
single dimensionless parameter*, B. In terms of the variables defined above, the buoyancy ratio 


can be written as 


Eel 
plu io a BPs bs 


pp & Pst Ps (Ps - Pst YH 
He ipo 


where Psp and Ps; are the air densities at the bottom and top boundaries of the computational 


domain, and H is the vertical dimension of this domain, as shown in figure 2. For typical low- 
sb - Pst 


Oo 


level atmospheric stratification, = 10> over a vertical range of H = O(100m). If the 


4 The steady motion of a plume with initial circular cross section in a uniform wind, U, depends on four 
independent dimensional variables: R, g, pp, and ae. Only one dimensionless parameter can be formed from 
v4 


Pp /R 


these four variables is, B = ——.. 
dO ps/ Oz 
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plume size is of R = O(10m), and with deficient density be = 10°, the corresponding buoyancy 
Oo 


ratio is B = 100. For large plumes, the expected range of the buoyancy ratio is J0 < B< 1000. 


We ran the computational code for four values of the buoyancy ratio, B = °, 50, 25, and 


digas pee aa neutral to strong stratification, with all other conditions: the initial plume size 
(Ry, Rz) = ( ), the initial deficient plume density ae = -].0, the initial height H7/R = 23.0, 


Va Va = 
and the buoyancy Reynolds number R,, = 1000, being the same. Symmetry across the plume 
centerline was enforced at all times. The discretization lengths for the computational elements is 
h/R =0.025 for the plume-air interface, and h,/R =0.5 in the stratified atmosphere, as shown in 
figure 2. The marching step in the wind direction is At =0.025°5. The vertical extent of the 
domain is 15<z <40, so that H/R=25. The horizontal extent is O< y <7.09. The 
computational results were used, besides investigating the fundamental dynamics governing 
plume rise and dispersion, to determine the values of the entrainment coefficient and the added 
mass coefficient in the trajectory formula, equation (7). It should be noted that the values of 
these constants may vary with B. Similar to Part I, we define the plume height at any downwind 
location as the algebraic mean of the maximum and minimum heights of the plume cross section. 

The computed plume center height of each case was fitted to the trajectory formula by 
adjusting the arbitrary constants, B and ky. The computed trajectories of the four cases, and the 
corresponding fitted curves are shown in figure 3. The best fit for all four cases were obtained | 
using the following pair of constants: B = 0.49 and k, = 0.7. It was found that while B had a 
stronger influence on the equilibrium height, k, affected the downwind distance where the 
equilibrium height was reached. As mentioned Benict equation (7) is expected to be valid only 
for the part of the trajectory prior to reaching the level-off height. 

The top curve in figure 3 corresponds to the case of plume rising in a neutrally stratified 


atmosphere. In this case, the plume rises continuously and there is no equilibrium height. 


5 In terms of the parameters used to check the convergence of the computations, such as the plume 


trajectory, numerical results with smaller values of h, hs and At are indistinguishable from the results 
described in this paper. 
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Equation (7) describes this curve using N2 = 0, B = 0.7, and k,= 0, as shown in Part I. 
Therefore, the same trajectory can be represented by the two-thirds law using two different sets 
of constants; one with the a finite added mass coefficient and the other without®. 

The lowest curve in figure 3, where the plume reaches the equilibrium height early, 
corresponds to a case with strong stratification. Weak, fast decaying oscillations are observed in 
the level-off region. This is consistent with the experimental observations that the plume 
oscillations decay rapidly in the downwind direction (Briggs, 1975). Plume oscillations are 
generated as the buoyant plume overshoots its equilibrium height and disturbs the initially calm 
background. These oscillations, which possess a well-defined frequency, then travel through the 
atmosphere in the form of gravity waves due to the restoring action of the stratified background. 
A detailed discussion on the origin and the character of these oscillations is provided in the 
Appendix using a simple analytical model. The intermediate cases also follow the extended two- 
thirds law well until the plumes reach their corresponding equilibrium heights, where the 
analytical curves overestimate the heights. As expected, the trajectory formula (4) generally 
cannot predict the maximum plume rise. 

Part of the kinetic energy associated with the plume motion, transformed from its initial 
potential energy through the vorticity, is radiated away in the crosswind section via the internal 
waves. Thus, a significant part of the surrounding air may be disturbed and a larger 
computational domain may be needed to capture the effect of these waves. To find out the effect 


of the finite domain size on the results, we repeated the calculation of the B= 25 case after 


6 Briggs (1969) proposed using B=0.4 in the formula for the plume radius, equation (5), and B=0.6 in the trajectory 
formula, equation (4). This means that, in the analytical model used to derive the two-thirds law, R = Ro + 0.6 z is 
used to evaluate the plume momentum, while R = Ro + 0.4 z , is used to evaluate the plume size. Using larger 
values of R in calculating the momentum may account for part of the momentum of the background fluid which is 
set in motion by the plume. Thus, it is not surprising that the early version of the two-thirds law, which did not 
incorporate the concept of added mass, could still describe the plume trajectory well (Fay et al, 1970, Hewett et al, 
1971). 
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doubling the horizontal extent of the computational domain to y =14.18. No significant changes 
were observed in the result up to x = 17.5 where the calculations were terminated. 

Now we compare the entrainment relation expressed by Equation (5) with the simulation 
results. In dimensionless form, if we use the equivalent radius’ to describe the size of the plume 


cross section, equation (5) can be written as 
Req = IN + Bz (9) 


Figure 4 shows the computed equivalent radius as a function of the plume rise. Clearly, the 
growth of the plume cross section depends weakly on the buoyancy ratio before reaching the 
level-off height. Beyond z=2, and before the level-off distance, the computed equivalent radius 
grows almost linearly with height. Both equation (9) and its corresponding form for a point- 
source plume, Reg = Bz, are shown using B = 0.49 which was obtained from the trajectory. As . 
shown in the figure, equation (9) overestimates the equivalent radius, although the slope is 
reasonably captured. Beyond z=2, the equivalent radius increases with height as if the plume 
had started from a virtual point source at z=0. 

It is interesting to note how the buoyancy ratio affects the plume spread close to the level- 
off height. In the neutrally stratified case, B = , the equivalent radius and the plume rise 
increase indefinitely. In the weakly stratified case, B= 50, the radius continues to expand after 
the plume reaches the level-off region. In the cases with stronger stratification, B= 25 and 12.5, 
the plume radius and rise are curtailed beyond the equilibrium height, as shown by the crowd of 


points near the level-off height. 
4.2. Other integral characteristics of the plume 


To investigate the effect of stratification on plume dispersion, we plot the evolution of the 
vertical and horizontal widths of the plume cross section, both as defined in Part I. Figure 5 
shows that beyond the level-off height, the vertical width of the plume cross section increases in 


7 The equivalent radius is defined as the radius of the circle whose area the same as the actual plume cross- 
sectional area. 
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the neutrally stratified case, oscillates around a nearly constant value in the weakly stratified 
atmosphere, and decreases in the strongly stratified atmosphere. Achieving an almost constant 
vertical width under conditions of intermediate stratification is consistent with experimental 
results (Briggs, 1975). The decrease of the vertical width under strong stratification is because 
while the top side of the cross section is capped, the bottom side continues to rise, as will be seen . 
in more detail in the next section. On the other hand, as shown in figure 6, the plume horizontal 
width increases continuously, with the rate of expansion being approximately the same for all 
cases. Thus, atmospheric stratification affects the plume dispersion mainly by suppressing the 
vertical fluid motion. 

As mentioned earlier, the effect of atmospheric stratification can be understood in terms 
of the vorticity generated in the background, as presented in figure 1. During the early stages, 
negative vorticity (on the right side of the centerline), which lifts the plume upwards, is 
generated at the interface between the plume and the surrounding. In the later stages, after the 
atmosphere has been disturbed, positive vorticity (on the same side), which tends to push the 
plume downwards, is generated in the stratified background. Using the results of the simulation, 
we can probe into this mechanism in more detail. Figures 7 and 8 show the total circulation, 
integrated over half of the computational domain, in the plume-air interface and in the 
background, respectively. Initially, and in all cases, the magnitude of the total circulation at the 
plume-air interface increases rapidly at a rate which is independent of B. However, The terminal 
value of the total circulation in the plume cross section is a function of B. The saturation of the 
total plume circulation is due to the reason described in Part I: both negative and positive 
circulation grows at the same rate following the large eddy roll-up stage. On the other hand, 
after the atmosphere has been disturbed, the positive background circulation increases 
continuously. The stronger the stratification, the faster the rate at which ambient circulation is 
generated. With strong stratification, the magnitude of the positive background circulation could 


be larger than that of the negative plume circulation. However, the induced motion by the former 
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is not necessarily stronger because this positive background circulation is distributed over an area 
much larger than the plume cross section, as will be shown in the next section. 

The deformation of the plume cross section can also be measured in terms of the growth 
of the boundary (interface) length between the plume and the surrounding atmosphere. The 
plume circumference, shown in figure 9, exhibits exponential growth following a short 
acceleration period with a high rate during the early large eddy roll-up stage and a lower rate 
during the small eddy roll-up stage. It is interesting to note that this circumference grows 
continuously beyond the point where the plume has reached a level-off height. This means that 
the buoyancy-generated turbulence still plays an important role in dispersing the plume material 
even at this later stage. The figure shows that the growth of the circumference is almost 
independent of the background stratification. 

The plume equilibrium, or level-off, height, Zeg, is an important quantity. This height is 
commonly used in conventional models as the terminal rise of a buoyant plume in the presence 
of atmospheric stratification. At this height, the plume material is assumed to be distributed 
according to a Gaussian distribution. Figure 10 shows the plume terminal rise as a function of 
the buoyancy ratio (two additional cases, B = 16.67 and 33.3, were computed to establish this 
curve). In practice, it is customary to use equation (4) to calculate the maximum plume height 
and multiply the result by an empirical factor to obtain the final equilibrium height. This, 
however, may not always be appropriate due to the following reasons: (a) the plume trajectory 
formula, equation (4), does not predict an accurate maximum plume height since the entrainment 
and similarity assumptions become invalid near this region; (b) the empirical multiplication 
factor is not a well defined quantity; and (c) the experimental results on point source buoyant 


plumes still involves fairly large uncertainties. 
4.3. Concentration and vorticity distributions 


Figure 11 shows the deficient density distribution across the plume cross section at 


successive downwind distances for the case of B = 12.5. The increment between the density 
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isolines is the same in each cross section and decreases for larger downwind distance. The 
number of isolines in the background increases for larger x because, while the number of isolines 
are the same in each section, the range of density variation decreases. Initially the plume cross 
section is circular (only half of the section is shown) and the background stratification is 
horizontally uniform, so that the background density isolines are all horizontal. Due to 
buoyancy, a double-vortex structure, which induces the upwards motion and deforms the plume 
cross section into a kidney-shaped pattern, is generated. Concomitantly, the initially horizontal 
background density isolines are disturbed and positive background vorticity (on the same half of 
the plume) is generated. 

Beyond x=7, each side of the cross section splits into two lumps with more material 
being associated with the lower lump. Both lumps continue to rise until the plume center reaches 
its equilibrium height at x=/2. Thereafter, the top lump stays almost at a fixed height while the 
lower one continues to rise slowly thus decreasing the overall vertical width of the plume. This 
bifurcation phenomenon was also observed in the two-dimensional thermal experiment of Tsang 
(1971). 


Figures 12 and 13 show contours of: (a) the plume deficient density, which is 


as ; (b) the vorticity; (c) the background deficient 
(4) 


proportional to plume material concentration, 


Ps 


zp ; and, (d) the velocity field in the cross-wind section at x = 4 and x = 1/4, 
(4) 


air density, 


respectively, all for the strongly stratified case with B= 12.5. Since we have used separate 
transport elements for the plume perturbation density and air density, these two fields can be 
calculated and presented separately here. 

At x = 4, the plume cross section exhibits the kidney shaped structure which was also 
observed in the neutrally stratified case (Part I). At this section, located only a short distance 
away from the source, negative vorticity is generated over a large portion of the plume-air 
interface, as demonstrated in figure 1. The cumulative action of this vorticity is to propel the 
plume upwards. Meanwhile, positive vorticity is generated in the surrounding, although its 


effect on the plume trajectory is not yet noticeable, see figure 3. The mechanism of positive 
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vorticity generation in the surrounding atmosphere can be explained using figure 12c which 
exhibits the density isolines in the immediate vicinity of the plume. Heavy air from lower 
altitudes is entrained into the core, and is displaced by lighter air which flows downwards from 
higher altitudes along the plume sides. This imparts strong curvature on the background density 
isolines and leads to the generation of positive vorticity. The velocity field, shown in figure 12d, 
depicts a structure which can be attributed to a vortex dipole, coincident with the plume material 
concentration contours. At x = 14, and according to figure 3, the plume is at the level-off region. 
On each side of the centerline, the plume breaks up into two lumps with the lower one containing 
more of the original plume material. The break-up of a plume into several lumps has been 
observed by Richards (1963) and is considered to be one of the important mechanisms of 
turbulent plume dispersion. 

The vorticity field in figure 13b shows two eddies: the lower eddy having strong negative 
circulation and the upper eddy having weak positive circulation (The sense of rotation of these 
eddies is confirmed by figure 13d). The positive background vorticity is distributed over a much 
larger area than that occupied by these two plume eddies. Since the plume is not rising at this 
stage, the net effect of these eddies on its trajectory must be balanced by the effect of the positive 
vorticity generated in the background. According to the conventional argument, the plume 
reaches the level-off height when its average density, being that of the original plume material 
after getting diluted by the entrained air, is the same as the local atmospheric density. This 
"dynamic" argument, however, relies on quantities, such as the plume size, which are difficult to 
define. Our "kinematic" argument is based solely on the motion induced by the vorticity 
generated within the plume cross section and in the background. 

The background density field in figure 13c illustrates the effect of the plume break-up on 
the motion of surrounding air. In the early stages, air, which is heavier than the plume material, 
is entrained vertically upwards into the core. At x=/4, the entrained air is forced to move 


horizontally away from the centerline, indicating that the plume has stopped its vertical ascent 
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and is spreading out horizontally. Meanwhile, some of the heavier air in the plume core is 
leaking downwards. 

The velocity field, depicted in figure 13d, depicts the manifestation of the complex 
vorticity structure within and in the immediate vicinity of the plume. As mentioned earlier, the 
vorticity field is composed of a positive eddy on the upper side and a negative eddy on the lower 
side of the nlyntig cross section. Near the center, y=0, the flow induced by this vorticity 
distribution acts to dislodge some of the air entrained earlier. It is important to add at this point 
that even after reaching the level-off height, the plume continues to disperse horizontally due to 


the action of the large scale vortices. 
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5. Conclusion 


The computational plume model, which is based on the utilization of the transport 
element method to solve the equations governing a buoyant flow in a horizontal uniform wind, 
has been extended to include the effect of atmospheric density stratification. Results have been 
obtained for a wide range of the buoyancy ratio which characterizes the relative importance of 
the plume buoyancy to the change in the atmospheric density. The computed plume trajectory 
and dispersion patterns have been compared with experimental measurements whenever 
possible. The extended two-thirds power law is found to describe the plume trajectory 
reasonably well before the plume reaches its equilibrium height. The entrainment and the added 
mass coefficients are found to be B = 0.49 and ky = 0.7, over a wide range of buoyancy ratio. 

The computational results suggest that the plume cross section first changes from a 
circular into a kidney-shaped form, and then breaks up into several lumps at later stage. This 
behavior has been confirmed qualitatively by various field and laboratory observations. 
However, quantitative measurements of the plume material distribution, which can be used in 
direct comparison with the numerical results, are currently scarce. This problem may be 
overcome using the newly developed rapid-scanning Lidar technique which is capable of 
obtaining the instantaneous density distribution within the plume cross section (Bennett et al, 
1992). Application of this technique may provide considerable field data to guide further model 
validation. 

Many features of buoyant plume motion, e.g., the kidney-shaped cross section, the 
induced motion of the surrounding air and the inhibiting effect of atmospheric stratification, can 
be attributed to the buoyancy-generated vorticity. This turbulence, being neither homogeneous, 
nor isotropic or stationary, occurs often in fluid motion with finite density gradient. This 
vorticity exhibits itself in the form of a long-preserved, large-scale, coherent vortex structure, 
which governs the entrainment process. Since this turbulence is frequently encountered in 


nature, its understanding and proper description is important both fundamentally and practically. 
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The agreement between the simulated results and the experimental observations shows 
the potential for applying the model to more complicated problems, such as, plume interaction 
with an inversion layer, heavy particle separation from buoyant plumes, and chemically reactive 
plume dispersion. Neglected from the current model are the initial momentum and/or shear in 


the plume section. Ways of modeling their effects are currently being investigated 
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Appendix 
On the level-off region of plume trajectory 

In a stable atmosphere, the air density decreases with height. As a buoyant plume rises, it 
entrains ambient air and the plume volume-averaged density increases. When the plume reaches 
the maximum height, it loses its vertical momentum since most of its kinetic energy is 
transformed into gravitational potential energy. In the level-off region of a plume trajectory, a 
_ slight overshoot of the final height is often seen. However, more than one oscillation is 
infrequently observed (Briggs, 1975). The same strong damping of oscillations is exhibited by 
our numerical simulations. To explain this behavior, Warren (1960) studied the damping of 
oscillation of a rising axisymmetric thermal and suggested that the rapid damping of plume 
oscillation is due to “wave drag”. In the following, we develop a simple model to explain this 
phenomenon. 

We assume that the plume motion is steady; the ambient wind, U, is uniform and the 
atmosphere is linearly stratified. The plume is assumed to be a parcel with volume-averaged 
density and at the position of its mass center. We denote the coordinate in the wind direction by 
x*, and the vertical distance above the equilibrium height by z*. The vertical momentum 
equation, after incorporating the Boussinesq approximation, can be written as 

pean Pe. (Al) 

oO x* Po 
where w* is the plume vertical velocity. The rate of change of the average density inside the 
plume is approximated as 

(pa Pp) =" (A2) 
where p, and pp, are the air density and plume deficient density, respectively; us is the influx of 
air per unit cross-sectional area into the plume due to entrainment. For a plume with larger 


deficient density, the entrainment is stronger, so that the influx of air into its core is higher. We 


parameterize this influx using an entrainment coefficient a: 


Fi ie OPp 3 (A3) 


26 


Realizing i <<1 and that the air density is vertically stratified, we can simplify equation (A2) 


a 


and rewrite it as 


0 Pp x 2 Pa _ 
Vacca ern lt, (A4) 


Eliminating p, from (A1) and (A4), we obtain 


owe Owe 
U2 + ou + N2 we =0 (AS) 
0 x*2 Ox* 
where N2=--& 2P5 is the square of the Brunt-Vaisaila frequency which characterizes the 


Po dz* 


stability of the atmosphere. Writing the vertical displacement of plume from the equilibrium 


height as E° = | w* dr*, then (A5) becomes 


2,* 
Tateinn ase te NCaetie 6 (A6) 


The solution of (A4) gives the trajectory of the plume: 


E (at) = bp exp (- Apa) cos Hh af 1-285 x). (A7) 


where the origin of x* is chosen at the location where E° reaches its maximum value, E, ; 
Clearly, this maximum amplitude is related to the vertical momentum of the plume when it first 
crosses the equilibrium height. 

Discussion: 

(1) The oscillation amplitude decays exponentially as the plume proceeds in the wind direction. 


The decay length scale is 


LE, =4U (A8) 
a 

which is longer for large ambient wind speed U or small entrainment coefficient a. Lg depends 

on N since @ is a function of stratification although it does not appear explicitly in (A8). For 


large N, or strong vertical stratification, the atmosphere is stable, and we expect the plume 
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entrainment to be weak, @ to be small and L, to be long. Therefore, the decay length scale is 


longer for stronger ambient stratification. 


(2) The spatial wavelength is 


N horas | (A9) 


which is longer for faster ambient wind speed, shorter for stronger stratification and weaker 
entrainment. 
(3) The number of plume oscillations depends on the ratio of 3Lg, the distance over which 95% 


of the maximum amplitude decays, and 27Ly, the wavelength of the oscillation, 


__3lq + 3u2Neal __a2 Al 
‘i 21d we el MD: are 


For longer decay length or shorter wavelength, more oscillations are expected. There are two 
extremes. (a) @=0, i.e., no entrainment is involved, the plume moves as a pure, non-dissipative 
internal wave. In this case, which corresponds to a passive plume in an internal wave field, we 
have an infinite number of oscillations. (b) @ =2N, very large entrainment is involved and no 


oscillations occur. It is interesting to note that n is independent of the wind speed. 
Generally, for 0 < a/2N < 1, the number of oscillations is finite and small if entrainment 


is large. The dependence of n on a@/2N is shown on figure Al. 
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Figure Captions 


Figure 1. Schematic drawing showing the generation of vorticity along the plume interface and 
in the surrounding, when atmosphere stratification is finite. 


Figure 2. The computational domain in the cross-wind section at time t = 0, showing the 
discretization of the plume-air interface and the surrounding into transport elements. The air 
density profile is shown on the right-hand side. 


Figure 3. The trajectories of the plume center in a linearly stratified atmosphere. Open symbols 
depict the numerical simulation results for four different buoyancy ratios. The extended two- 
thirds law is shown by solid lines for the same cases. 

Figure 4. The equivalent radius of the plume cross section versus the plume rise for different 
buoyancy ratios. Numerical results are shown by open symbols while the integral relations are 


depicted by solid lines. 


Figure 5. The change of the plume vertical width, defined in terms of the maximum vertical 
extension of the plume cross section, along the downwind distance. 


Figure 6. The variation of the plume horizontal width, described by the maximum crosswind 


horizontal extension, along the downwind distance. 


Figure 7. The total circulation generated at the plume-air interface on the right half of the plume 


cross section as a function of the downwind distance. 


Figure 8. The total circulation generated in the background on the right half of the plume cross 
section as a function of the downwind distance. 


Figure 9. The length of the plume circumference along the downwind distance. 


Figure 10. The terminal or equilibrium height of a buoyant plume rising in a linearly stratified 


atmosphere as a function of the buoyancy ratio. 
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Figure 11. The evolution of the plume cross section, depicted in terms of the density contours, at 
different downwind locations for the case of B=12.5. All frames are from y=0.0 to y=7.0 and 
2=23.0 to z=30.0. 


Figure 12. Contours of (a) the plume density, (b) vorticity, (c) the background density and (d) 
the velocity vectors in the plume cross section at x = 4 for a strongly stratified atmosphere, 
B=12.5. All frames are from y=0.0 to y=6.0 and z=24 to z=30. 

Figure 13. Contours of (a) the plume density, (b) vorticity, (c) the background density and (d) the 
velocity vectors in the plume cross section at x = 14 for a strongly stratified atmosphere, B=12.5. 


All frames are from y=0.0 to y=6.0 and z=24 to z=30. 


Figure Al. Number of oscillations near the equilibrium height as a function of a/2N . 
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